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Summary 

The LEWICE ice accretion computer code has been extended to include the solution of the two- 
dimensional Navier-Stokes equations. The code is modular and contains separate stand-alone pro- 
gram elements that create a grid, calculate the flow field parameters, calculate the droplet trajecto- 
ry paths, determine the amount of ice growth, calculate aeroperformance changes, and plot 
results. The new elements of the code are described. Calculated results are compared to experi- 
ment for several cases, including both ice shape and drag rise. 


Nomenclature 


A + = 26 

A p = Droplet surface area 
c p = Specific heat 
Cj = Drag coefficient 
Cj = 5 
C 2 = 2000 
D - Drag 

E = Convective flux terms in component 
of momentum equation 
E m = Total collection efficiency 
f = Freezing fraction 
F = Convective flux terms in rj- compo- 
nent of momentum equation 
g = Gravitational constant 
h c = Convective heat transfer coefficient 
i - enthalpy 
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J = Jacobian of coordinate transforma- 

tion 

k s = Roughness height 

Lf = Heat of freezing 

L v = Heat of vaporization 

l = Mixing length 

m = Mass 

Q = Solution matrix 

q = Heat transfer 

r c = Recovery factor 

S = Viscous flux terms in Navier-Stokes 

equation 

s = Surface distance 

T = Temperature 

u z = Sheer velocity 

V = Velocity, Volume 

x = Spatial coordinate 

y = Spatial coordinate 

a = Angle of attack 

(3 = Local collection efficiency 

k = Clauser constant, 0.0168 


l 



|i, = Eddy viscosity 

v = Kinematic viscdsity 

p = Density 

T] = Transverse coordinate in body-fitted 

system 

£ = Stream-wise coordinate in body fitted 

system 

Subscripts 

a = Air 

c = Convection, Incoming from cloud 
e = Evaporation, Edge of boundary layer 

I. Introduction 

Over the past several years the need for ice accretion prediction capabilities has been growing. 
Many aircraft and ice protection system manufacturers have used the NASA developed ice accre- 
tion code, LEWICE. 1 The LEWICE code predicts the growth of ice on 2D surfaces through appli- 
cation of an inviscid panel method module, a particle trajectory calculation module, and a control 
volume energy balance/ice growth calculation module. The use of this code assists the analyst in 
assessment of potential icing hazards due to ice growth on unprotected surfaces and in the appro- 
priate placement of candidate ice protection systems. Several upgrades to this code are currently 
underway to improve the capability of LEWICE. Some of these include; extension to 3D geome- 
tries, 2 inclusion of a thermal ice protection system model, 3 and the addition of a performance deg- 
radation evaluation capability. 4 

An additional area of potential improvement for the LEWICE code comes from the introduction 
of an alternate flow code calculation method. The current LEWICE code uses the inviscid panel 
method, S24Y, developed by Hess and Smith. 5 This approach has the advantage of being a very 
fast code and produces good results over an adequate range of conditions. There are several prob- 
lems associated with the use of this code which have been documented in the LEWICE User’s 
Manual 1 and by other users of the code.^ The S24Y code uses a Karmen-Tsien compressibility 
correction. This helps extend the range of applicability to higher Mach numbers. However, the ba- 
sis for this correction is a small perturbation approximation, hence thick airfoils and high angles 
of attack may not be accounted for using the correction. Additionally, the convective heat transfer 
values, which contribute significantly to the energy balance/ice growth calculation, may be quite 
different under conditions of high Mach number or of large local pressure gradients due to surface 
curvature. 

Another area of concern is the evaluation of velocities near the surface at panel edges. The calcu- 
lated velocities at the panel edges can become much larger than actual values, resulting in the pos- 
sibility of unrealistic reversals in water droplet trajectories. This problem has been addressed in 
the LEWICE code by the creation of an artificial surface located some distance upstream of the 
actual surface. The droplets, now avoiding the flow reversals, impinge on this pseudo-surface and 
the local collection efficiency is calculated based on this surface. 

The collection efficiency value is a measure of the amount of water impinging on the surface com- 
pared to the amount of water passing through a plane upstream of the airfoil bounded by the air- 
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foil thickness. The local collection efficiency is a measure of how droplet impingement varies 
along the surface of the airfoil. The local collection efficiency is based on an evaluation of the 
droplet impingement pattern as a function of distance along the surface. Thus, the pseudo-surface 
approximation can markedly change the collection efficiency values if it differs significantly from 
the actual surface contours. This has the potential for problems when determining the collection 
efficiency on iced surfaces. 

In order to include more accurate velocity results near the airfoil surface and to include the effects 
of compressibility, the inviscid panel method was replaced by an Euler code. The use of an Euler 
code in place of the panel code in LEWICE was documented in a previous article. This change 
also requires the use of a grid generator and the particle trajectory calculation must be altered to 
use interpolated velocities from a grid as opposed to determining the velocities by evaluating the 
contribution of each panel to the potential function at the desired location. The grid generator used 
has the ability to develop body-fitted grids for extreme shapes, such as glaze ice horns. The veloc- 
ity interpolation calculation is incorporated into a subroutine within the particle trajectory code 
module. 

Due to the irregular shapes associated with ice formations, multiple stagnation points can develop 
on the iced surface. This can lead to difficulty in selecting the point at which to start the boundary 
layer calculations for the upper and lower surfaces. LEWICE deals with this situation by either 
asking the user to select one of the stagnation points to start the boundary layer calculations or to 
create another pseudo- surface encompassing the stagnation points and calculating a single new 
stagnation point for that surface. The use of a grid based code can lead to different approaches to 
the solution of this problem. For example, the velocities at grid point locations away from the sur- 
face could be used to determine which grid line divides the upper surface from the lower surface 
and the surface point associated with that grid line could be the starting point for the boundary 
layer calculations. Such an approach could be automated and eliminate the need for user interac- 
tion as is currently the case in LEWICE. In this investigation, the point closest to the stagnation 
point from the previous time step is used to start the integral boundary layer calculation. 

II. Code Description 

The code consists of four major modules tied together by a simple command-line user interface. 
The four code modules are a grid generator, an Euler/Navier-Stokes flow solver, a particle trajec- 
tory calculation, and an energy balance/ice growth calculation. The code modules produce graphs 
of pressure coefficient distributions, lift and drag coefficient histories, residual histories, particle 
paths, collection efficiency distributions, flow field properties at the edge of the boundary layer, 
thermal conditions on the surface and the calculated ice shape. Contour plots can be produced 
which show the conditions in the far field. Output listings from each of the modules are also pro- 
duced in order to allow examination of exact numerical values for parameters of interest. Files for 
geometry and flow solution information are written in PLOT3D** format in order to allow compat- 
ibility with other grid generators and flow solvers which could be substituted for those used in this 
code. 

II-l. Grid Generation Module 

The grid generator used in LEWICE/NS is a hyperbolic equation solver developed by Barth. ^ 
This code creates a C-grid around the user defined surface and can be easily applied to complex 
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surfaces such as an iced airfoil. The number of grid points and grid point spacing along the sur- 
face can be adjusted to concentrate grid points in regions of high curvature. This gives the user 
considerable flexibility in trying to develop an acceptable grid. An example of the grid generator 
output is shown in Figure 1, the grid for an iced NACA0012 airfoil. 

The equations solved in this grid generator are hyperbolic and take the form; 


x^+y^ = ° Hq. 1 

H y r\~ X r\y^ = J = V Eq>2 


where x and y are the Cartesian coordinates and £, and rj are the coordinates of the body-fitted sys- 
tem. The inverse of the Jacobian, J, approximates the local cell volume, V. Since this is a set of hy- 
perbolic equations, they are solved using a marching procedure starting at the body surface and 
moving along the q direction out to the outer boundary. 

One reason for choosing this grid generator is that it has been modified to work for surfaces with 
concave curvature or surface slope discontinuities. This is especially helpful when trying to re- 
grid the geometry after the ice shape has been added. This grid generator has been used on several 
complex shapes as documented by Barth/ 1 It has also been used extensively^’ 1 ^ for ice accretion 
geometries and has produced usable grids, for most cases with minimal user interaction. The issue 
of grid spacing in conjunction with the accuracy of the flow solution can be a problem even for 
clean airfoil geometries, however careful monitoring of the flow solution for decrease of the resid- 
uals, convergence of integrated results such lift, and development of inappropriate flow results 
(e.g. excessive entropy generation) should allow a high degree of confidence in the solution. 

Since the overall LEWICE/NS code is modular, an alternate grid generator can be easily substitut- 
ed for the current module. The only criterion for substitution is that the grid generator output the 
grid coordinates in PLOT3D format. This allows the grid file to be used by the other modules as 
well as by the PLOT3D visualization package. 

II-2. Euler/Navier-Stokes Flow Solution Module 

The Euler/Navier-Stokes code used in LEWICE/NS is the ARC2D code developed by Steger iq 
and Pulliam. 11 This code uses the grid file and an input file which describes the flow conditions 
and selects some code procedure parameters. The code solves the Euler or Navier-Stokes equa- 
tions in the body-fitted coordinate system. The equations solved are, 

d x Q + + d^F = Re-'d^S Eq. 3 

where the details of the terms in the equations can be found in references 10 and 11. The Cartesian 
coordinates ( x,y ) and the body-fitted coordinates (^,T|) are obtained from the grid generator. De- 
tails of the numerical algorithm used to solve these equations are also found in references 10 and 
11 . 

In addition to the flow solver, a turbulence model must also be used when the code is run in the 
Navier-Stokes mode. The turbulence model currently used in ARC2D for icing calculations is an 
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algebraic eddy viscosity model developed by Potapczuk.^ This model does not require the deter- 
mination of the boundary layer edge as in the Cebeci-Smith model nor does it require the determi- 
nation of values along grid lines as in the Baldwin-Lomax model. This model is based on the 
mixing length concept and defines the eddy viscosity through appropriate length and velocity 
scales. In this model, the eddy viscosity, p,, is given by the equation, 

p, = p/ 2 |co| Eq- 4 


In this model, the length scale is determined by the conditions at the surface and then saturates at 
some distance from the surface as indicated in the data of Klebanoff. 13 Hence, the length scale is 
given by the relation. 


and, 
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Eq. 5 


Eq.6 


where the constants Cj and C 2 are determined empirically by matching the p, and l contours, for 
flat plate turbulent boundary layers, from the calculations with experimental results (normally C] 
= 5 and C 2 = 2000). The above equations include the terms y* and y + which are measures of the 
transverse dimension within the boundary layer and are defined as. 


y + = y/y 


where. 



Eq. 7 


Eq. 8 


and where t is the shear stress at the wall. 

The mixing length and vorticity can be determined for any grid point with no need to march out 
from the surface along grid lines. 

In addition to the turbulence model, a method of modeling roughness must be added to this calcu- 
lation. The rough surfaces of an ice shape can cause premature transition from laminar to turbu- 
lent flow. This in turn can alter the drag and heat transfer characteristics of the iced airfoil. The 
Cebeci-Chang 14 model was chosen for use with ARC2D along with the turbulence model de- 
scribed above. The Cebeci-Chang model represents the roughness as a displacement to the bound- 
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ary layer that is factored into the calculation of the mixing length. The amount of the 
displacement. Ay, is based on the roughness height k 5 and is calculated from the relation 


Ay + = 


0.9 


^“<ex p 


f -kZ 


5 < kZ < 70 


0.7 (O 


0.58 


70 < k+<, 2000 


Eq. 9 


where Ay + = Ayu x /v and k s + - k s u x /v. 

The mixing length is then found by adding Ay into Eqs. 5 and 6. This produces the following 
equations. 


y + Ay 
* 
y 


<c x 


l (y) = 


Ci * 

K c: y 


/ 

( ,y + Ay ^ 

V * / 

1 


1 - 

i- y 


(l-e y 

c i 


V J 

V 

V J 

J 



Eq. 10 


and, 


y + Ay 
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>C, 


l (y) = K-^y* 

c 2 


Eq. 11 


The velocity and pressure values resulting from the flow field calculation are used for several pur- 
poses in the LEWICE/NS code. They contribute to the evaluation of the particle trajectories. The 
values at grid lines away from the surface are used in the determination of boundary layer edge 
conditions. These edge values are used in the integral boundary layer calculation of the energy 
balance/ice growth module. Finally, the velocity and pressure field is used to determine the lift 
and drag values for the iced airfoil following the ice growth calculation. 

The user determines the validity of the flow field solution by examining residual and lift histories 
and the pressure coefficient distribution and entropy contours. The flow field solution values are 
contained in the Q matrix, which is stored in PLOT3D format for plotting and to insure a standard 
file transfer process. 

Normally, the flow solver is used in Euler mode during the ice accretion portion of the calculation. 
This allows the code to be run in a faster manner due mostly to the fewer number of grid points. 
There are some conditions for which the Navier-Stokes solution is required. This usually occurs 
when there are regions of high surface curvature resulting in large pressure gradients. These large 
gradients can result in convergence problems with the Euler equations. The viscous terms in the 
Navier-Stokes equations can introduce sufficient dissipation to allow convergence. 

Once the final ice shape has been determined, the flow solver can be run again with the viscous 
terms included. This calculation provides information on the changes in lift and drag as a result of 
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the ice shape. This usually requires the use of the grid generator again in order to develop a grid of 
sufficiently fine resolution to capture the flow gradients near the surface. 

II-3. Particle Trajectory Module 

Once the flow field for a given geometry has been obtained, the velocities are used to determine 
the trajectories of super-cooled water droplets from release points upstream of the body to impact 
locations on the surface. In the original LEWICE code, this was done with a Lagrangian predic- 
tor-corrector method in a code developed for LEWICE. 1 This code is also used for LEWICE/NS, 
however, some changes were required in order to deal with a grid-based flow code instead of the 
original inviscid panel method, which can determine velocity values at any arbitrary location. 

The equations of motion for water droplets are based on the assumption that they are rigid spheres 
and hence the only forces acting are drag and gravity. This assumption is considered valid for 
droplet radii of less than 500 pm. 1 The governing equations are thus. 


where 


m'x = - D cosy-t- mgsina 
my = - D siny— mg cos a 


Eq. 12 


y = 



Eq. 13 


D = 



Eq. 14 


V = Ju„-V x ) 2 +(y p -V t ) 2 


Eq. 15 


and where V x and V y are the components of the flow field velocity at the droplet location and 
x p , y p , x p andy^ are the components of the droplet velocity and acceleration, respectively. 

The water droplet trajectories are calculated using an Adams-Boulter predictor-corrector algo- 
rithm. After each calculation of a droplets location, it is checked to determine if it has impacted on 
the body. The calculation continues until a droplet impacts the surface or until it has passed some 
user selected location, indicating that the droplet has missed the body. Details of this calculation 
are found in reference 1. 

As mentioned in the introduction, the original LEWICE code requires a pseudo- surface in order to 
avoid unrealistic velocities near the surface which lead to inaccurate droplet trajectories. The use 
of a grid based code avoids this problem. Any inviscid flow solution does require some special 
treatment near the surface. The slip flow at the surface is not physically correct and thus the drop- 
let trajectory calculation does not use those values. Instead, once the particle has crossed the first 
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grid line off of the surface, the droplet path is no longer altered by the flow field and is considered 
to be tangent to the previously calculated path line. The intersection of that path line with the body 
surface is taken as the impingement location. 

The pattern of droplet impingement on the body surface determines the amount of water that hits 
the surface and becomes part of the ice growth process. The ratio of the actual mass that impinges 
on the surface to the maximum value that would occur if the droplets followed straight-line trajec- 
tories, is called the total collection efficiency, E m . The total collection efficiency for the body is 
found by integrating the local collection efficiency, P, between the upper and lower limits of drop- 
let impingement. The local collection efficiency is defined as the ratio, for a given mass of water, 
of the area of impingement to the area through which the water passes at some distance upstream 
of the airfoil. Taking a unit width as one dimension of both area terms, the local collection effi- 
ciency can then be defined as. 


n _ dy o _ A y 0 
” ds As 


Eq. 16 


where Ay 0 is the spacing between water droplets at the release plane and Ay is the distance along 
the body surface between the impact locations of the same two droplets. The local collection effi- 
ciency is illustrated in Figure 2. 

The local collection efficiency is the necessary input for the energy balance/ice growth module. It, 
along with the free stream velocity and the cloud liquid water content, determines how much wa- 
ter impinges on the local region of the surface under consideration. Variations in the local collec- 
tion efficiency can significantly alter the ice growth for that surface region. 

n-4. Energy Balance/Ice Growth Module 

The growth of ice on the surface is a complex fluid dynamics, heat transfer, and mass transfer pro- 
cess. The incoming water may freeze on impact or some fraction may freeze while the remaining 
water either runs along the surface or collects in pools. The processes determining which of these 
occur include surface tension effects, roughness, skin friction between water and ice or water and 
airfoil surface, shear forces between water and air, and convection and conduction heat transfer. A 
simplified model of this process has been developed by Messinger 1 *’ and is used in the original 
LEWICE code. This model is used in the LEWICE/NS model as well. However, as alternate ice 
growth models are developed, they may easily be substituted into this code due to its modular na- 
ture. 

The Messinger model, as implemented in the LEWICE code, is described fully in Appendix A of 
the LEWICE User’s Manual. 1 As such, only a brief description will be given here. The ice growth 
process is modeled as a control volume analysis. The control volume is bounded by the body sur- 
face and by an arbitrary boundary considered to be at the edge of the boundary layer. The two 
chordwise boundaries coincide with constant ^ grid lines established by the grid generation code. 
The lower boundary of the control volume is initially the body surface but it moves outward with 
the ice growth process, remaining at the top of the solid-fluid interface. For dimensional com- 
pleteness, the control volume is considered to extend one unit length in the spanwise direction. 
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The control volume is used to perform a mass and energy balance, as depicted in Figures 3 and 4. 
The equation governing the mass balance is. 


m r + m r —m. — m. = m ; 

C 'in C ' out 1 


Eq. 17 


where m c is mass flow rate of incoming water, m r is the mass flow rate of runback water from 
the previous control volume, m e is the mass flow rate of evaporated water, m r is the mass flow 
rate of water running back to the next control volume, and m i is the mass flow°rate of water leav- 
ing the control volume due to freeze-out. The equation governing the control volume energy bal- 
ance is. 


m c l w,T +m rJw, sur(i- 1) Eq jg 

= Wv. sur + "V Jw, sur + "Vi, sur + Vc As + ^k As ) 

where z' w> T is the stagnation enthalpy of the incoming water droplets, i w sur ^_ j) is the enthalpy 
of the water flowing into the control volume from upstream, i v sur is the enthalpy of the vapor 
leaving the control volume due to evaporation, i w sur is the enthalpy of the water running back to 
the next control volume, z ( - sur is the enthalpy of the ice leaving the control volume, q c is the heat 
transfer due to convection, and q ^ is the heat transfer due to conduction at the bottom of the con- 
trol volume. 

The incoming energy due to water droplet impingement and runback are calculated from known 
information. The energy leaving the control volume due to evaporation and convection can be cal- 
culated independently. Currently, the convection heat transfer calculation is performed in the 
same manner as in the original LEWTCE code. That is, an integral boundary layer analysis is per- 
formed using the Euler calculation results as the boundary layer edge conditions. The code could 
be used in Navier-Stokes mode to obtain these values, but the integral boundary layer calculation 
is faster and thus was used in this investigation. The heat transfer due to conduction is not consid- 
ered in this analysis, as the ice layer is considered to act as an insulating surface. This leaves the 
energy loss due to freeze-out and the energy leaving the control volume due to runback as un- 
knowns. In particular, the mass flow rates for these two terms are unknown, as was the case for 
the mass balance. This leaves two equations and two unknowns and the system can be solved. The 
details of the evaluation of each of the terms in the energy equation can be found in Appendix A 
of the LEWICE User’s Manual. 1 

A useful concept for evaluation of the nature of the ice accretion being calculated is the freezing 
fraction. This is the fraction of the total water coming into the control volume that changes phase 
to ice. The equation defining freezing fraction is, 


f = 


m : 


+ m. 


Eq. 19 


This term can also be used to simplify the evaluation of the energy balance. 
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The convection heat transfer term plays an important role in the LEWICE/NS energy balance. It is 
through this term that the aerodynamics and the roughness levels can influence the development 
of the ice accretion. Currently, the convection heat transfer is determined from an evaluation of 
the boundary layer growth on the surface, using an integral boundary layer method. The pressure 
distribution determined by the Euler calculation is used as input to the boundary layer calculation. 
The boundary layer calculation determines the displacement thickness and the momentum thick- 
ness. The Reynold’s analogy is used to determine the heat transfer coefficient. Roughness is ac- 
counted for by a correlation developed by Ruff. 1 The complete description of the integral 
boundary layer calculation is found in Appendix B of the LEWICE User’s Manual. 1 

Expanding the terms in the energy equation as described in the LEWICE manual and combining 
Eqs. 17-19 yields the following form of the energy equation, 

vi~ 

*c[^. <7,- 273.15) + T 

+,A rJ c f ..,.,„_„< T ,u,(i- 1 ) - 273 . 15 )] +q k &s 
= m e [c p _JT„ r -m.l5) +LJ 

Eq. 20 

+ [0-/> 0*e + ">r h ) -™e^ Pw ,JT sur -213.\5) 

+/(m c -m r J[c / , iw (7’, Hr -273.15)-L / ] 

rV 2 ~ 

+ h c T sur -T e -j^ As 

Pa _ 

The solution process for this equation is started by identifying the stagnation point. Since no run- 
back water can enter the control volumes on either side of the stagnation point, one of the un- 
knowns is determined for these two control volumes. Thus, the solution may be marched back, on 
the upper and lower surfaces, towards the trailing edge from the stagnation point. As m i is deter- 
mined for each control volume, the resulting ice growth thickness can be found from the ice den- 
sity and the dimensions of the control volume. The ice thickness values define a new iced airfoil 
geometry by adding that thickness to the body in a direction normal to the surface. In regions of 
high curvature, the new ice surfaces can intersect or diverge.The method for dealing with the re- 
definition of the surface under these circumstances is the same as that of the original LEWICE 
code and is described in Reference 1 . Once the new geometry has been defined, the entire process 
can be started again at the grid generation step. 

III. Comparison to measured ice accretion and drag 

The LEWICE/NS code was used to calculate ice shapes on a NACA 0012 airfoil for a series of 
conditions tested by Olsen, et al, 16 in the NASA Lewis Icing Research Tunnel (IRT). In that test 
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program, the ice shapes and resulting drag values were measured for a series of tests on the 
NACA 0012 airfoil with varying icing conditions. The code was used to try and reproduce the re- 
sults of the test for both ice shape and drag. 

m-1. Ice shape predictions 

The Olsen tests covered a wide range of conditions and examined such effects as temperature, 
drop size, and liquid water content (LWC). The most dramatic effect uncovered, with respect to 
drag, was the correlation of drag with temperature. Figure 5 shows the variation of drag as a func- 
tion of temperature. At temperatures just below freezing, there is a large increase in drag which 
tapers off as the temperature decreases. The suspected reason for this is the change from a glaze 
ice condition to a rime ice condition. The glaze shapes, with their associated horns, cause a signif- 
icant change to the pressure distribution over the airfoil resulting in the large change in drag. 
Rime ice, on the other hand, has a smaller effect on the pressure distribution with the drag rise be- 
ing mostly due to skin friction. 

This series of sixteen temperature conditions (with LWC x V x Time held constant) is detailed in 
Table 1. The LEWICE/NS code was used to evaluate these cases in order to evaluate its ability to 
predict ice accretion over a wide range of conditions. There were some difficulties in obtaining 
solutions for some of these cases. These problems were either with grid creation or with conver- 
gence of the flow field calculation. Further investigation is required in order to determine what 
steps need to be taken to overcome these difficulties. In all, eleven of the fourteen cases with ice 
depositions were calculated. Measured coordinates for two of these eleven were not available and 
so only nine cases were compared. The results of the LEWICE/NS analysis indicates reasonable 
agreement with the rime ice conditions and fair to poor agreement with the glaze ice conditions. 
Additionally, the transition from rime to glaze shapes appears to occur at different temperatures in 
the code results than in the experimental results. This can be seen in Figs. 6 and 7. 

The differences between calculation and experiment can be attributed to two sources. The 
LEWICE/NS pressure distributions and droplet trajectories are not the same as those of the origi- 
nal LEWICE code. This was demonstrated in a previous investigation. 7 As a result, the distribu- 
tion of incoming water on the surface is different between the two codes as is the convective heat 
transfer acting on that water. The differences in these two calculations leads to significant differ- 
ences in the resulting ice shapes. The roughness correlation used in LEWICE, which has a signif- 
icant effect on the convective heat transfer coefficient, was tuned to the LEWICE results in order 
to produce reasonably accurate ice shapes over a limited range of icing conditions. The results of 
this investigation indicate that the correlation is not appropriate for the LEWICE/NS code. This 
means that either a new correlation is required for the LEWICE/NS code or the convection heat 
transfer must be obtained from the energy equation directly. The former approach would require a 
parametric investigation over several icing cloud conditions, such as temperature, velocity, LWC, 
or droplet size, to determine correlations between input roughness levels and output ice shapes. In 
this approach, the development of a suitable correlation is dependent on how representative the 
correlating data set is of all possible icing conditions. However, it does have the advantage of pro- 
ducing satisfactory results over some range of conditions. The latter approach would necessitate 
the use of the viscous option in the LEWICE/NS code during each time step of the ice accretion 
calculation. 
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Without considering the correlation between prediction and experiment, there are some interest- 
ing features of the LEWICE/NS ice shape results which deserve further comment. In one case, 
shown in Fig. 8, the ice shape produces a shadow region for droplet impact. That is, a protuber- 
ance in the ice shape blocks water droplets from impacting on a portion of the surface aft of the 
protuberance. Droplets originating from points lower along the release plane do impact further aft 
on the surface, once they pass under the protuberance. This results in two distinct collection effi- 
ciency curves as seen in Fig. 9. 

The convective heat transfer coefficient in this region is such that the freezing fraction is less than 
unity. This means that some of the water on the protuberance travels back into the region with no 
impingement. The water then rapidly freezes before traveling back to the region where impinge- 
ment resumes. In the aft impingement region, the heat transfer is high enough to cause rapid 
freezing and a large horn develops in this region as seen in Fig. 6(c). This secondary horn is also 
evident in the experimental shape corresponding to this condition, albeit of smaller size. This in- 
terplay of collection efficiency and convective heat transfer is postulated to be a possible mecha- 
nism for the formation of ice fingers, seen frequently on actual ice shapes but never actually 
modeled. A similar horn is also seen in the Fig. 6(b), again for both the calculated and experimen- 
tal results. 

III-2. Iced airfoil drag predictions 

Results from previous investigations by Shin, et al, 4 have indicated very good agreement between 
calculations and experimental values. Their calculations captured the large rise in drag at warm 
temperatures associated with glaze ice conditions. There were some differences in the level of 
drag which they associated with transition and roughness modeling. In a later investigation, 17 
they were able to account for this difference by extending the rough region well beyond the ice 
growth region on the airfoil. In an effort to avoid using this artifice, the LEWICE/NS code is be- 
ing used to examine the drag results of the calculated ice shapes described in the previous section. 

Since many of the calculated ice shapes did not match the experimental values, especially for the 
glaze ice conditions, it was not expected that the resulting drag values would compare favorably. 
However, the drag values were evaluated in order to determine if the large drag rise at warm tem- 
peratures would be seen in the calculations. The drag value for these calculations were produced 
by using the standard LEWICE roughness value as the value for k s used in Eq.9 (i.e. k s = 
XKINIT). The k s values used in the calculations are shown in Table 1. These roughness levels 
were applied only over the region with ice on the surface. A value of k s = 0.5x10' 5 was used in the 
clean region of the surface. This was done in order to simulate some minimal roughness level of a 
clean surface. The resulting drag values are shown in Fig. 5. 

Drag calculations for glaze ice geometries have been calculated before using this roughness mod- 
el. 12 The results have indicated that the pressure drag is a much more influential term than the vis- 
cous drag for glaze ice formations with substantial horn growths. As a result, it is expected that 
the drag predictions for rime ice formations will be more sensitive to the roughness and turbu- 
lence modeling described above. 

The drag results parallel the ice shape results in terms of agreement with experimental values. The 
low speed results agree quite well with the measured drag values at low temperatures. At warmer 
temperatures, where the pressure component of the drag is predominant, the results do not com- 
pare as well. The agreement between calculated and measured ice shapes is also not that good at 
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these conditions. The high speed results do not agree as well as the low speed results. The drag is 
underpredicted at all temperatures. The drag values at warmer temperatures were not obtained for 
this velocity because of difficulties in determination of the final ice shape, as described previously. 
It can be seen, though, that even at -12°C the drag values begin to differ markedly between the 
calculation and experiment. Comparison of the ice shapes for this condition, Fig. 7(b), indicates 
that the experimental condition is a glaze formation while the calculated result is a rime shape. 
This explains the marked difference in drag results for this case. As a result of the difficulties in 
determining the ice shape for the warmer conditions, the ability of the code to reproduce the char- 
acterisitic drag rise could not be evaluated. 

IV. Concluding Remarks 

A new version of the LEWICE family of codes has been developed which substitutes the use of a 
Navier-Stokes code for the potential flow code, S24Y. The code is normally used in an Euler 
mode during the ice accretion portion of the calculation with the viscous terms added only during 
the follow-on performance degradation calculations. As documented previously, 7 the new code 
has the capability to calculate ice shapes for high Mach numbers and large angles of attack. 

Comparisons of code results to experimental ice tracings were good for rime ice conditions but 
were only fair to poor for glaze ice accretions. The differences are considered to arise from poor 
correlation of the roughness parameter with the boundary layer edge velocity and pressure distri- 
butions calculated by the Euler portion of the code. Close correlation between the LEWICE code 
roughness parameter and the convective heat transfer coefficient was necessary to produce good 
agreement between ice accretion prediction results and actual ice tracings in the original 
LEWICE. This correlation will have to be adjusted to allow use of the Euler code within the 
LEWICE framework. Development of an alternate correlation suggests that some additional ex- 
perimental data on heat transfer, collection efficiency, and roughness characterization for ice 
shapes would be useful for enhancing the robustness of any such correlation. An alternative 
would be to utilize the solution of the energy equation, obtained from a Navier-Stokes calculation 
prior to each ice accretion calculation time step, in place of the integral formulation currently 
found in LEWICE and LEWICE/NS. This latter approach may still require correlation informa- 
tion, however it may be on a lower level of approximation, such as the heat transfer coefficient 
over a characteristic roughness element. Further investigation of the relationship between rough- 
ness level and ice shape thus seems warranted no matter which approach is pursued. 
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Table 1: Test Conditions for Icing Tests with NACA 0012 Airfoil, Chord = 

0.53m, Angle of Attack = 4° 
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Figure 1 - Near field grid for NACA 0012 airfoil with ice 
shape on leading edge. 



Figure 2 - Definition of local collection efficiency. 
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Figure 3 - Mass balance control vol- 
ume. 


Figure 4 - Energy balance control vol 
ume. 
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Figure 5 - Comparison of measured and calculated drag 
coefficent results for several temperature conditions 
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Figure 6 - Measured and calculated ice shapes for V=209 km/hr, LWC=1.3 g/m , 
MVD=20jJ.m, Duration = 8min (except calculated (d) = 6min). 
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(c) T tot = -8°C 


■ Clean Airfoil I I Measured Ice Shape I I Calculated Ice Shape 

Figure 7 - Measured and calculated ice shapes for V=338 km/hr, LWC=1.05 g/m 3 , 
MVD=20jim, Duration = 6.2min. 
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Figure 8 - Droplet trajectories from LEWICE/NS showing shadow 
region between main ice growth and horn. 
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Figure 9 - Collection efficiency distribution resulting from 
droplet trajectories of Figure 8. 
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